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Abstract. 

We describe an axisymmetric general relativistic code for rotational core collapse. The code evolves the coupled 
system of metric and fluid equations using the ADM 3 + 1 formalism and a conformally flat metric approximation 
of the Einstein equations. Within this approximation the ADM 3 + 1 equations reduce to a set of five coupled non- 
linear elliptic equations for the metric components. The equations are discretized on a 2D grid in spherical polar 
coordinates and are solved by means of a Newton-Raphson iteration using a block elimination scheme to solve the 
diagonally dominant, sparse linear system arising within each iteration step. The relativistic hydrodynamics equa- 
tions are formulated as a first-order flux-conservative hyperbolic system and are integrated using high-resolution 
shock-capturing schemes based on Riemann solvers. We assess the quality of the conformally flat metric approxi- 
mation for relativistic core collapse and present a comprehensive set of tests which the code successfully passed. 
The tests include relativistic shock tubes, the preservation of the rotation profile and of the equilibrium of rapidly 
and differentially rotating neutron stars (approximated as rotating polytropes), spherical relativistic core collapse, 
and the conservation of rest-mass and angular momentum in dynamic spacetimes. The application of the code 
to relativistic rotational core collapse, with emphasis on the gravitational waveform signature, is presented in an 
accompanying paper. 
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1. Introduction 

The advent of gravitational wave astronomy - major in- 
terferometric gravitational wave detectors have already 
started taking data - aggravates the need for gravitational 
wave templates computed from theoretical models of po- 



tential sources of gravitational radiation (Pradier et al 
2001). Among the most promising sources are aspherical 
core collapse supernovae. However, reliable core collapse 
simulations, in particular those of rotational core collapse, 
require a relativistic treatment of the dynamics because 
of the counteracting stabilizing effect of rotation and the 
destabilizing effect of relativistic gravity. In addition, the 
simulations should incorporate a realistic equation of state 
(EoS) and an accurate treatment of neutrino transport. 

Due to these complexities, all previous investigations 
have either neglected or approximated one or all of 
these requirements (see, e.g. Miiller (1998 ) and references 
therein). Up to now there exist no successful multidimen- 
sional numerical simulations of rotational core collapse to 
neutron stars in general relativity, even for simple mat- 
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ter models. On the other hand, simulations which include 
better microphysics in the form of realistic (nuclear) equa- 
tions of state or neutrino transport have either been con- 
fined to spherical symmetry or have used Newtonian grav- 
ity. 

There is an obvious need to improve this situation 
if one wants to answer questions such as: What is the 
quantitative influence of general relativity on the dynam- 
ics of core collapse compared to Newtonian gravity? How 
do higher densities and velocities expected in relativistic 
simulations change the properties of the rotating proto- 
neutron star? To what extent do relativistic effects alter 
the gravitational wave signal? How does the rotation rate 
evolve in a relativistic simulation, and what consequences 
does that have on the development of dynamical or secular 
instabilities in the neutron star? 

To this end we have started a project to compute re- 
liable gravitational wave templates from rotational core 
collapse and nonspherical core collapse supernovae. In a 
first step we have improved the treatment of the dynam- 
ics beyond the still commonly used Newtonian approach 
and computed the general relativistic hydrodynamics of 
rotational core collapse consistently with the evolution of 



2 



H. Dimmelmeier et al.: Relativistic rotational core collapse. I 



the axisymmetric dynamic spacetime. The computational 
tool developed by us for this purpose is described in this 
publication, while its application to the gravitational col- 
lapse of rotating stellar cores approximated as polytropes 



is devoted to an accompanying paper ( Dimmelmeier et al 
2002()7]We note that in pimmelmcicr et al. (2001| ) we have 



already presented preliminary results of our investigation. 
General relativistic simulations of improved models in- 
cluding both a realistic equation of state and an approx- 
imate treatment of weak interactions and neutrino trans- 
port are planned in the near future. 

The asumptions we adopt in the present and accom- 
panying paper are as follows: Our matter model obeys an 
ideal gas EoS with the pressure P consisting of a poly- 



with a being the lapse function, which describes the rate 
of advance of time along a timelike unit vector nor- 
mal to a hypersurface, (3 l being the spacelike shift three- 
vector which describes the motion of coordinates within a 
surface, and jij being the spatial three-metric. 

In the 3 + 1 formalism, the Einstein equations split into 
evolution equations for the three-metric 7y and the extrin- 
sic curvature K y , and constraint equations which must be 
fulfilled at every spacelike hypersurface: 



-2aK lJ +W t (3 j +W j (3 l , 



d t Kij = -ViVj-a + a(Ri 



KKij - 2K lk Kf 



tropic and a thermal part (Janka et al. 1993; Zwerger & 



Jij 



Miiller 1997). Since we are mostly interested in the gravi- 
tational radiation emission, which is controlled by the bulk 
motion of the fluid, we neglect any other microphysics 
like electron capture and neutrino transport. The initial 
configurations are rotating relativistic polytropes in equi- 
librium ( jKomatsu et al. 1989a| Jb"[ |Stergioulas fe Friedman 



R + K 2 - K i3 K lj - 167T/9H, 
Vi(K ij - j ij K ) - 8irS j . 



(2) 



(3) 

(4) 
(5) 



In these equations Vj is the covariant derivative with re- 
spect to the three-metric 7y, is the extrinsic curva- 
ture, Rij is the Ricci tensor, R is the scalar curvature and 
1995|)~}vhich are marginally stable. Their collapse is ini- K — K\ is the trace of the extrinsic curvature. The matter 



tiated by decreasing the effective adiabatic index 7 P to 
some prescribed fixed value. Our initial data solver allows 
us to construct both uniformly and differentially rotating 
polytropes. 



fields appearing in the above equations, £y, S J and pu, 
are the spatial components of the stress-energy tensor, the 
momenta and the total energy, respectively. 

Following Wilson et al. (1996) (see also Flanagan 



We solve the general relativistic hydrodynamic equa- 
tions numerically. The metric equations are approximated 
by adopting the so-called conformal flatness condition for 
the three- metric 7y ( Wilson et al. 1996| ). This results in a 
constrained system of elliptic metric equations which are 
much simpler than the full system of Einstein equations. 
The dimensionality of the problem is reduced by assum- 
ing both axisymmetry with respect to the rotation axis, 
and equatorial plane symmetry. The numerical code uses 
a spherical grid in spherical polar (r, 9) coordinates. 

This paper is organized as follows: In Section we 
present the mathematical framework used in our ap- 
proach, namely the conformally flat metric equations and 
the relativistic hydrodynamic equations. Section ^ is de- 
voted to describing the initial models. Next, in Section [| 
we describe all relevant issues concerning our numerical 
implementation. Section [| presents a comprehensive num- 
ber of tests we have used to validate our code. The paper 
ends with a summary in Section ^. 

Unless otherwise stated we use geometrized units (G = 
c = 1). Greek (Latin) indices run from to 3 (1 to 3). Latin 
indices are raised and lowered using the three-metric. 



2. Mathematical framework 

2.1. Metric equations 



(1999[ ); [Mathews fc Wilson (2000| )) we approximate the 



general metric g M „ by replacing its spatial three-metric 
jij with the conformally flat (CF) three-metric (confor- 
mal flatness condition - CFC hereafter): 



Jij 



where 7^ is the flat metric (7;. 



(6) 

in Cartesian coor- 



dinates). In general, the conformal factor </> depends on 
the coordinates (t,r,8). Therefore, at all times during a 
numerical simulation we assume that all off-diagonal com- 
ponents of the three-metric are zero, and the diagonal el- 
ements have the common factor <fi 4 . Note that all metric 
quantities with a hat are defined with respect to the flat 
three- metric 7y. 

Within this approximation the ADM equations reduce 
to a set of five coupled elliptic (Poisson-like) equations for 
the metric components, 



A0 - -2ir^ [ phW 2 -P+ ^-^1 

\ 167T / 



(7) 



We adopt the ADM 3 + 1 formalism ( [Arnowitt et al. 1962| ) 
to foliate the spacetime into a set of non-intersecting 
spacelike hypersurfaces. The line element reads 



/ 7K K l i\ 
A(a</>) = 2™</> 5 ( ph(3W 2 - 2) + 5P + , (8) 

\ 167T J 

A/3* = leira^S 1 + 2K i ^ j (^j - ^V l V k p k , (9) 

where V and A are the flat space Nabla and Laplace op- 
erator, respectively. The transformation behavior between 
the extrinsic curvature defined on 7^ and 7^ is as follows: 

-*Ku, K ij = <l>- w K ij . 



K, 



ds 2 



2 dt 2 + ^{dx 1 + (3 i dt)(dx j + (3 j dt), 



(1) 



,j Y " — V -11 . (10) 

The metric equations (^-^|) couple to each other via 
their right hand sides, and in case of the equations for fi % 
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via the operator A acting on the vector j3 l . The equa- 
tions are dominated by the source terms involving the 
hydrodynamic quantities p, P, and v l , whereas the non- 
linear coupling through the other, purely metric, source 
terms becomes only important for strong gravity. On each 
time slice the metric is solely determined by the instanta- 
neous hydrodynamic state, i.e. the distribution of matter 
in space. 

Applying the CFC to the traceless part of Eq. (g) yields 
the following relation between the conformal factor and 
the shift vector components: 



6 



(11) 



With this, Eq. (Q) transforms into a definition for the ex- 
trinsic curvature components: 



(12) 



The above approximation is exact in spherical symme- 
try. Nevertheless we have used it for our simulations of 
rotational core collapse. The rationale behind this choice 
is the well-justified assumption that the energy emitted 
in form of gravitational waves is only about 10 -6 of the 
energy released during gravitational collapse in a Type II 



supernova event (Muller 1982; Monchmeyer et al. 1991 
Zwerger fc Muller 1997| ). Furthermore, during core col- 



lapse deviations from spherical symmetry are modest even 
for rapidly rotating cor es compared to tho se occuring in 
a neutron star merger ( IWilson et al. 1996| )). This fact is 
also reflected in the strength of the gravitational wave sig- 
nal from core collapse, which is three orders of magnitude 
smaller than that resultin g fro m the plunge phase of a 
compact binary. In Section 5.5 we will discuss the quality 
of the CFC in the context of rotational core collapse in 
more detail. 

As a consequence of setting the off-diagonal ele- 
ments of jij to zero, the degrees of freedom represent- 
ing gravitational waves are removed from the spacetime. 
Therefore, gravitational wave emission is calculated in 
a post-processing step using the quadrupole formula. A 
more detailed description of the wave extraction methods 



is given in the accompanying paper (Dimmelmeier et al 
2002j)~| 

2.2. General relativistic hydrodynamic equations 

The hydrodynamic evolution of a relativistic perfect fluid 
with rest-mass current J M = pu^ and energy-momentum 
tensor T^ v = phu^u" +Pg tJ ' 1 ' in a (dynamic) spacetime g^ v 
is determined by a system of local conservation equations, 
which read: 



V M J" = 



0. (13) 

Here is the covariant derivative and it M is the four- 
velocity of the fluid. Its three-velocity, as measured by an 
Eulcrian observer at rest in a spacelike hypersurface £t is 



+ 



Following the work of Banyuls et al. (1997) we now in- 
troduce the following set of conserved variables in terms of 
the primitive (physical) hydrodynamic variables (p, Vi,e), 
with p and e being the rest-mass density and specific in- 
ternal energy density, respectively: 



D 

T 



phW 2 v\ 



J» n = p hW 2 



D. 



(15) 
(16) 
(17) 



In the above expressions, is the unitary vector nor- 
mal to the slice and \- l v denotes the projection operator. 
Furthermore W is the Lorentz factor defined as W = cm 1 



and satisfying the relation W = l/\/l — ViV 1 . In addition, 
h = 1 + e + P/p is the specific enthalpy and P is the 
pressure. 

The local conservation laws, Eq. (|l3|) , can be written as 
a first-order, flux-conservative hyperbolic system of equa- 
tions, 



1 



-9 



dx° 



dx" 1 



= Q 



(18) 



with the state vector, flux vector and source vector given 

by 



u 



Q = 



[D,S J: r] 
D [ v l - 



, Sj 



a 



5)P, 



v — 
0, T^ v 



£ I + /V 

a 



9g, 



d In a 

dx^ 



- 1 u V 9\j 



(19) 



(20) 



(21) 



Here y/—g — a^/j, with g — det(g M „) and 7 — det(jij) 
being the determinant of the four- metric and three-metric, 
respectively. In addition, J 1 * are the Christoffel symbols. 



3. Initial models for rotational core collapse 

3.1. Equation of state 

For a rotating iron core before collapse the polytropic re- 
lation between the pressure P and the rest mass density 
P, 



P = Kp 1 



(22) 



with 7 = 7 ini = 4/3 and K - 4.897 x 10 14 (in cgs units) 
is a fair approximation of the density and pressure strati- 



fication (Zwerger & Muller 1997) 



To model the physical processes which lead to the onset 
of collapse, we reduce the effective adiabatic index from 
7ini to 71 on the initial time slice. During the infall phase 
of core collapse, we also assume the matter of the core to 
obey a polytropic EoS, Eq. (|22|), which is consistent with 
the ideal gas EoS for a compressible inviscid fluid, 



(14) P=( 7 _l) /0e , 



(23) 
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where e is the specific internal energy. 

During and after core bounce, when the shock starts 
to propagate out, the matter in the post-shock region is 
heated, i.e. kinetic energy is dissipated into internal en- 
ergy. Therefore, in this stage we demand that the to- 
tal pressure consists of a poly tropic and a thermal part, 
P = P p + Pth- We further assume that the polytropic in- 
dex jumps from ji to 72 for densities larger than nuclear 
matter density p nuc . This approximates the stiffening of 
the EoS. Requiring that the pressure and internal energy 
are continuous at the transition density p n uc one finds for 
the polytropic internal energy 



,2/3-^^2 



K 



,71-1 



72 — 1 



n 72-l „7l — 12 
r rnuc 



(72 ~ li)K 
(72 - l)(7i - 1) 



71-1 



(24) 



for p < pnuc and p > p nuc , respectively (see Zwcrgcr 
(1995|)]- The thermal part of the pressure, which is due 
to shock heating, is given by Pth = (7th — l)/° e th, where 
7th = 1.5 and e t h = e — e p - This describes a mixture of 
relativistic (7 = 4/3) and nonrelativistic (7 = 5/3) com- 
ponents of an ideal fluid. 

From these considerations one can construct an EoS 
for which both the total pressure P and the individual 
contributions P p and P t h are continuous at p nuc : 



P = 



7 ~7th 
7-1 

+ (7th - 



Kpl^P 1 



(7th ~ 1) (7-7i) 
(71 - 1)(72 - 1) 



Kp 



nuc P 



l)pe. 



(25) 



3.2. Rotating relativistic stars in equilibrium 

Our initial models are obtained using Hachisu's self- 



consistent field method (Komatsu et al. 1989a) which is 
a general relativistic method for solving the hydrostatic 
equilibrium equations for rotating polytropic matter dis- 
tributions. The general metric to describe rotating ax- 
isymmetric relativistic matter configurations in equilib- 
rium is 



3 2 = -e 2p dt 2 + e 2& (dr 2 + r 2 d9 2 ) 

„2/§„2 -2 Cl{ Ai *, , ,A4\1 



+e 2/ V sin 2 6 {dip - udt) 2 , 



(26) 



with the metric potentials v, a, and u). 

The only nonvanishing three- velocity component is v 3 , 
the velocity of the fluid as measured by a ZAMO (zero 
angular momentum observer). Thus, the Lorentz factor is 
given by W — — V3V 3 . With the definitions u° = 

We"", and u 3 — QWe~" ' , where fl is the angular velocity 
of the fluid as measured from infinity, the four-velocity 
reads: 

= We'" {1,0,0, ty. (27) 
Thus, one gets for the cp-component of the three-velocity: 



e- p {VL-uj), 



r 2 sin 2 9{fl -co), 



V3 



v v = v V3V 3 = e 13 "r sin#(f2 — uj) 



(29) 
(30) 



Here v v is the proper rotation velocity with respect to 
the ZAMO. The specific angular momentum of the fluid 
j = u U3 is given by 

W 2 e 20~v) r 2 sin 2 e{Q _ u) = V_3*f (31) 



J 



(i-v 3 v 3 )(Q-uy 



Using the Einstein equations one can derive the equa- 
tion of hydrostatic equilibrium for an axisymmetric mass 
distribution rotating with the angular velocity D, = f2(r): 



VP 

ph 



1 — vi 



jwn = 0. 



(32) 



For barotropes where p depends only on p the integra- 
bility of Eq. (^|) requires that j is a function of Q, only. 
Hence, the simplest choice for a rotation law is 

j=A 2 {Q c -n), (33) 

where f2 r is the value of f2 at the coordinate center and 



A is a positive constant ( Komatsu et al. 1989a ). In the 
Newtonian limit the rotation law (|33|) reduces to 



n 



A 2 



1 



A 2 + d 2 



for A 
for A 







(34) 



= < A^_ 

I 2 

where d = rsmO is the distance from the rotation axis. 
The upper case in Eq. (34) corresponds to a rigid rotator 
and the lower one to a configuration with constant specific 
angular momentum. 

Using the rotation law specified in Eq. ( p3| ) one can 
integrate the equation of hydrostatic equilibrium ( |32| ) to 
obtain (Komatsu et al. 1989a) 



ln/i- 



\ ln(l - v%) - X -A 2 {^ - fl c ) 2 = const, 



(35) 



The four-metric components d, 7 = /3 + i>, 5 = v — and u> 
are given by the four coupled partial differential equations 



(Komatsu et al. 1989a) 

A (^/ 2 ) = S S , 



A 



A 

2 



1 



7e 



t/ 2 = S-, 



(36) 
(37) 

(38) 
(39) 



where A is the Laplacian and p = cos 6. The source terms 
Sg, Sj, Su, and Sa depend in general on all four met- 
ric functions. In order to solve this system of equations 
the first three PDEs are converted into integral equations 
using Green functions which in turn are expanded into 
Legendre polynomials. In discretized form, these equations 
are iterated until convergence is obtained. During each it- 



(28) eration, the discretized equation (39) for a is integrated. 
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Table 1. Set of simulated collapse models: A is the degree 
of differential rotation, /3 ro tini is the initial rotation rate 
(i.e., the initial ratio of rotational energy and the absolute 
value of the gravitational binding energy), and 71 (72) 
is the adiabatic index at densities below (above) nuclear 
density. See text for further details. 



Collapse 


A 


/^rot ini 


7i 


72 


model 


[10 8 cm] 


[%) 






A1B1G1 


50.0 


0.25 


1.325 


2.5 


A1B2G1 


50.0 


0.5 


1.325 


2.5 


A1B3G1 


50.0 


0.9 


1.325 


2.5 


A1B3G2 


50.0 


0.9 


1.320 


2.5 


A1B3G3 


50.0 


0.9 


1.310 


2.5 


A1B3G5 


50.0 


0.9 


1.280 


2.5 


A2B4G1 


1.0 


1.8 


1.325 


2.5 


A3B1G1 


0.5 


0.25 


1.325 


2.5 


A3B2G1 


0.5 


0.5 


1.325 


2.5 


A3B2G2 


0.5 


0.5 


1.320 


2.5 


A3B2G4 soft 


0.5 


0.5 


1.300 


2.0 


A3B2G4 


0.5 


0.5 


1.300 


2.5 


A3B3G1 


0.5 


0.9 


1.325 


2.5 


A3B3G2 


0.5 


0.9 


1.320 


2.5 


A3B3G3 


0.5 


0.9 


1.310 


2.5 


A3B3G5 


0.5 


0.9 


1.280 


2.5 


A3B4G2 


0.5 


1.8 


1.320 


2.5 


A3B5G4 


0.5 


4.0 


1.300 


2.5 


A4B1G1 


0.1 


0.25 


1.325 


2.5 


A4B1G2 


0.1 


0.25 


1.320 


2.5 


A4B2G2 


0.1 


0.5 


1.320 


2.5 


A4B2G3 


0.1 


0.5 


1.310 


2.5 


A4B4G4 


0.1 


1.8 


1.300 


2.5 


A4B4G5 


0.1 


1.8 


1.280 


2.5 


A4B5G4 


0.1 


4.0 


1.300 


2.5 


A4B5G5 


0.1 


4.0 


1.280 


2.5 



3.3. Initial models and collapse parameters 

For a given adiabatic index 7 ro tini, the initial models are 
determined by three parameters: The central density p c 
the degree of differential rotation A (see Eq. ([53])), and 
the rotation rate /3 r otini 7 which is given by the ratio of 
rotational energy and the absolute value of the gravita- 
tional binding energy. The central density of an iron core 
is about 10 10 g cm -3 . We adopt this value for p c - ln i in all 
initial models. 

Due to the lack of consistent models of rotating iron 



cores from stellar evolution calculations, we follow Zwerger 
fe Miiller (1997) and treat A and /3 r otini as free parameters 



of the initial models. The only other free parameters are 
71, the adiabatic index at subnuclear densities (p < p n uc)> 
and 72, the adiabatic index at supranuclear densities (p > 
p nuc ). As in Zwerger & Miiller (1997), we set 72 = 2.5 and 
71 equal to 1.325, 1.320, 1.310, 1.300 or 1.280, respectively. 
In total, we have evolved 26 models (see Table Q). One of 
the models (A3B2G4 so f t ) has also been run with 72 = 2.0 
in order to test the influence of a softer supranuclear EoS 



Table 2. Atmosphere parameters used in our core collapse 
simulations. For the extremely differentially rotating mod- 
els A4 the threshold value / a tmthr is in the range 5 x 10 -4 
to 1 x 10- 4 . 



fatm thr 


fatm 


Patm 

[g cm" 1 ] 


EoS Vi 


10 -5 


10 -io 


1 





on the collapse dynamics. The nuclear density is set to 
= 2.0 x 10 14 g cm -3 in all simulations. 



3.4. Atmosphere treatment 

One common problem of many hydrodynamic codes is 
their inability to handle regions where the density is zero 
(or very small compared to the typical densities of in- 
terest). Such regions are encountered, for example, when 
simulating pulsating stars on an Eulerian grid, or when 
simulating rotating stars on a spherical grid. In the first 
case a very low density "atmosphere" is required outside 
the stellar surface into which the star can expand dur- 
ing its pulsations. In the second case the stellar surface 
is flattened due to centrifugal forces, i.e. part of the com- 
putational grid extends outside the flattened surface and 
must be filled with some low density atmosphere. Thus, 
when computing initial models of rotating cores we intro- 
duce a low density atmosphere outside the stellar matter 
distribution (see Table ||). An atmosphere is assumed in 
all zones where the density p is below a prescribed thresh- 
old value which is practically defined as a fraction / a tm thr 
of the central density of the initial model, i.e. the density 
of the atmosphere p a t m — /atm/Ocim- During a simulation 
the atmosphere is reset after each timestep. This ensures 
that it adapts to the time-dependent "shape" of the stellar 
surface. 

For the values specified in Table |^ and an initial central 
density p C i n i — 10 10 g cm" 3 the homogeneous atmosphere 
has a density p^tm — 1 g cm" 3 . For all practical purposes 
this can be considered as a low density atmosphere, which 
has no influence on the dynamics of core collapse, the only 
tradeoff being a slight loss of angular momentum to the 
atmosphere from regions of the star which are close to the 
boundary (see Section 5.4). The EoS in the atmosphere is 
assumed to be polytropic, too. 

4. Numerical implementation 

4.1. Computational grid 

We use spherical polar coordinates (t, r, 9, ip), and assume 
axial symmetry with respect to the rotation axis, i.e. none 
of the variables depends on the azimuthal coordinate ip. 
However, motions in the (^-direction are allowed, and the 
associated velocity V3 and momentum S3 are nonzero. 

We further assume symmetry with respect to the equa- 
tor at 9 = 7r/2. Therefore, 9 spans the interval between 
and 7r/2. The angular grid consisting of ng zones is equally 
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spaced, with 6j denoting the coordinates of the cell cen- 
ters in angular direction (1 < j < n$). The n r radial zones 
(zone centers located at ) are logarithmically spaced. 
Each computational cell {ri,6j) is bounded by interfaces 
and 6*j+i, respectively. The angu- 



at r,- 



6, 



^ir/ng, while the radial size is 



lar size of a zone is A9 
given by Ari = r i+ i — r i _i. The radial size of the three 
innermost cells is enlarged in order to ease the timestep re- 
striction; for i > 3, Arj+i/ Ari = const. Note that r i = 
is the origin, and r nr+ i = R defines the outer bound- 
ary of the grid. The rotation axis coincides with 9i = 0, 
and 9 + i = i"7r denotes the equatorial plane. We use 
n r = 200 radial zones and ng = 30 angular zones, which 
corresponds to an angular resolution of 3° . The radial size 
Ari of the innermost cells is about 500 m. Convergence 
tests show that this grid resolution is sufficient to resolve 
the important flow features of core collapse dynamics. 



4.2. Boundary conditions 

We impose symmetry conditions for both hydrodynamic 
and metric quantities at the center (r = 0) and the rota- 
tion axis (9 = 0), and at the equatorial plane (9 = n/2). 
Quantities are assumed to be either symmetric or antisym- 
metric across a boundary. We combine these two options 
in the following notation for a quantity qij : 



where 



± = 



+ symmetric, 

antisymmetric. 



(40) 
(41) 
(42) 



(43) 



The symmetry conditions for the scalar quantities p, e, 
4> and a are straightforward. Like their Newtonian coun- 
terparts p, e and $ they have to be continuous across all 
boundaries. Thus, the scalar quantities are assumed sym- 
metric (+) with respect to all boundaries. Contrary to 
Newtonian hydrodynamics, the symmetry conditions for 
the vector quantities are nontrivial, as they occur both in 
covariant and contravariant form. 

By defining "physical" velocities (the same can be done 
for the shift vector components 



v r 
vg 



ViV 1 



Vl 

- 2 r-\ 2 



b 2 v\ 



j 2 2 
9 TV , 



V3V 



~ 2 r 1 sin 1 0V3 



"r sin 9v 



(44) 
(45) 
(46) 



which are 

\vg\Av, 



bounded by the speed of light (0 < 
\v r \, \vg\, \v v \ < 1), we can specify the symmetry condi- 
tions from Newtonian hydrodynamics. Keeping in mind 
the symmetry behaviour of the geometrical terms in the 
three- metric components 7^ and 7^ , the symmetry condi- 



tions for the "physical velocity vector" and the shift vector 
read: 





v r vg 


v v 


/3 1 


I3 2 /3 3 


center 








+ + 


pole 


+ - 




+ 


- + 


equator 


+ - 


+ 


+ 


- + 



Here we demand that (i) all velocities have to vanish at 
the center (no mass flow across the center), that (ii) the 
meridional velocity vg is zero along the rotation axis and 
in the equatorial plane (no mass flow across these bound- 
aries), and that (iii) the azimuthal velocity v v is zero along 
the rotation axis. Hence, these entries are set to (— ), while 
all others, where the velocities are continuous, are set to 
(+). The shift vector components /3 l are treated according 
to the contravariant three- velocity vector v l , as the shift 
vector corresponds to a "coordinate velocity" . 

At the outer radial boundary, the hydrodynamic quan- 
tities are extrapolated from the interior to the boundary 
zones. However, this procedure cannot be used for the met- 
ric quantities, as they arc determined by elliptic equations 
which define a boundary value problem. We thus deter- 
mine their boundary conditions by matching the interior 
metric to an exterior spacetime metric. We assume that 
the exterior spacetime is given by the Schwarzschild solu- 
tion for a spherically symmetric vacuum spacetime. The 
quality of this approximation is sufficiently good for all 
practical purposes provided the spacetime at the surface 
of the star does not deviate too much from spherical sym- 
metry. This is true for iron core initial data, and also holds 
for rotational core collapse, as the outer boundary of the 
radial grid is fixed (at the equatorial radius of the ini- 
tial model) at about 1500 km and as the monopole term 
always dominates the gravitational field. 

The Schwarzschild solution in isotropic coordinates is 



1 



ds z = - 



2r 



^2 _J_ -^gra v 



dt 2 



+ 1 + 



2r 

M. 



grav 



2r 



(dr 2 +r 2 dfl 2 ). 



(47) 



For a vanishing shift vector this metric can be exactly 
matched to the conformally flat interior metric: 



1 + 



2r 



a 



(1-%-) 
(1 + %-) 



f3 l = 0. (48) 



As the gravitational mass M grav is conserved and as f3 l is 
always close to zero at the outer boundary, this matching 
to the isotropic Schwarzschild metric provides an accept- 
able outer boundary condition for the metric coefficients 

(f> 7 a and (5 l . 

We stress that these boundary conditions are only an 
approximation to the exact vacuum spacetime solution, 
since no rotating spacetime can be described by a con- 
formally flat metric. By the above matching we neglect 
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effects like frame dragging (/3 3 ^ 0), which would in gen- 
eral be present in the vacuum spacetime region outside 
the computational domain. 

4.3. The hydrodynamic solver 

The hydrodynamic solver of our code performs the numer- 
ical integration of system (|l8| ) using a so-called Godunov- 
type (high-resolution shock-capturing - HRSC) scheme. 
Such schemes are specifically designed to solve nonlinear 
hyperbolic systems of conservation laws (see, e.g., Toro 



(1997 ) for definitions). In a Godunov-type method the 



knowledge of the characteristic structure of the equations 
is essential to design a solution procedure based upon ei- 
ther exact or approximate Riemann solvers. These solvers 
compute, at every cell-interface of the numerical grid, the 
solution of local Riemann problems. Therefore, they au- 
tomatically guarantee the proper capturing of discontinu- 
ities which may arise naturally in the solution space of a 
nonlinear hyperbolic system. 

The time update of system ( |l8| ) from t n to t n+1 is per- 
formed according to the following conservative algorithm: 



U 



n+l 



At 

An 
At 

A9 



( F uj + l/2 - F i,j-l/2) 



(49) 



The numerical fluxes F andF are computed by means of 



Marquina's approximate flux- formula ( Donat & Marquina 
199G| ). The fluxes are obtained independently for each di- 



rection. The time update of the state- vector U is done si- 
multaneously in both spatial directions using a method of 
lines in combination with a second-order (in time) conser- 
vative Runge-Kutta scheme. Moreover, in order to set up 
a family of local Riemann problems at every cell-interface 



we use the piecewise parabolic method (PPM) of Colclla 



& Woodward (1984) for the reconstruction of cell interface 



values, which provides third-order accuracy in space. 

Time derivatives of all metric quantities (0, a, and 
P l ) are required to compute the Christoffel symbols, and 
hence the source term Q. As the equations for the met- 
ric in the CFC approximation do not provide explicit ex- 
pressions for the time derivatives of these quantities, we 
approximate them numerically by discretized derivatives 
using values from the two time slices S(n and E t n-i, e.g. 



dt 



At n - 



(50) 



Finally we note that after the hydrodynamics update 
we need to recover the primitive variables (p, Vi , e) from 
the conserved ones (D,Si,r). Since the relation between 
these two sets is only given implicitely, we have to resort 
to a Newton-Raphson iterative method. This procedure is 
discussed in Appendix |A|. 



4.4. The metric solver 

One of the standard methods to solve an elliptic linear 
system like Poisson's equation on a numerical grid is to 
discretize the equation. Let the vector of unknowns be u = 
(uij), and the vector of sources be / = (fi,j), where 
labels the position (r j , 9j ) on the grid. This discretization 
leads to a linear matrix equation with the right hand side 
being the source vector of the original equation: 



Au = f. 



(51) 



In general the matrix A is sparse, i.e. the number of 
nonzero elements is much smaller than the total number 
of elements. For solving sparse linear systems there ex- 
ist a variety of efficient numerical methods and software 
packages for various computer architectures. 

Although the CFC metric equations ([t]-^) do not form 
a linear system like Poisson's equation, they can still be 
reduced to a linear system. Our strategy is as follows: We 
first define a vector of unknowns, whose components are 
the five metric quantities: 



(52) 



Then the metric equations can be written as 

f(u) = 0, (53) 

with / denoting the vector of the metric equations for u. 

The explicit form of the metric equations ( |53| ) is given 
in Appendix |b[ We discretize these equations on the (r, 9)- 
grid, denoting the metric components u k (k — 1,...,5) 
at the cell-center (rt,9j) by u t - •. Their values at cell- 
interfaces (needed for the hydrodynamic solver) are ob- 
tained by interpolation. The radial and angular derivatives 
are approximated by central finite differences. A 9-point 
stencil is required to formulate the discretized equations. 
Hence, the system of metric equations, when discretized, 
gives rise to a nonlinear system of dimension 5 x n r x ng , 



/(«) = o, 



(54) 



with the vector of discretized equations / = fij = ff*j 
for the discretized unknowns u = ui, m — u™ m . For this 
system we have to find the roots u\ j . 

We use a multi-dimensional Newton-Raphson itera- 
tion method to solve the nonlinear system, i.e. the non- 
linear system of dimension 5 x n r x ng is reduced to a 
linear one of the same dimension for the Jacobi matrix 
of /. This linear matrix equation has to be solved once 
in each iteration step. For our set of equations the Jacobi 
matrix of the nonlinear system is in general both sparse 
and diagonally dominated, which reduces the complexity 
of the linear problem significantly. 

The performance of our axisymmetric general relativis- 
tic code is mostly determined by the computation time 
spent in the metric solver, and thus by the computational 
efficiency of the linear solver. Therefore, the linear prob- 
lem has to be solved as efficiently as possible. After exten- 
sive investigation we have decided to use a block tridiago- 



nal sweeping method (Potter 1973). As the matrix which 
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15.0 



Fig. 1. Decrease of the maximum increment of the 
unknown vector Am^ 3X with the number of Newton- 
Raphson iterations s for a typical metric computation. 
Until the tolerance measure Au^ ax reaches the limit set by 
the machine precision at s i=s 8, it decreases exponentially. 
After that it saturates well below the precision threshold 
used in the code, which is marked by the horizontal dotted 
line. 

defines the linear problem has nine bands of blocks rather 
than three (due to the 9-point stencil; see above), we com- 
bine the left three, the middle three, and the right three 
bands, collecting them in three bands of submatrices - left, 
diagonal and right. Although these submatrices are sparse, 
the block tridiagonal sweeping algorithm leads to dense 
submatrices which have to be inverted. This is achieved 
using a standard LU decomposition scheme for dense ma- 
trices as implemented in the LAPACK numerical library. 

As far as computational performance is concerned, the 
block tridiagonal sweeping method exhibits an n r x rig 
behavior, and is thus superior to both the Gauss-Jordan 
elimination and the bi-conjugate gradient stabilized meth- 
ods which we have also tried. In particular, as it is a direct 
solver it does not suffer from convergence problems. 

The tolerance measure we use to control convergence of 
the Newton-Raphson method is the maximum increment 
of the solution vector on the whole computational grid: 



max (Am s ) = max (Au^J). 



(55) 



We are able to meet this convergence criterion for a toler- 
ance level of Au max < 10~ 15 in all our computations. 

In Fig. [I] we show how Au max decreases with the num- 
ber of iterations s in the metric computation of a typical 
core collapse model around the time of maximum density. 
The convergence measures exhibits an inverse exponential 
behavior until saturation is reached. This occurs when the 
increment of u s approaches machine precision. 

In addition, we monitor the residual of the function 
vector ry = f(u s )) on the entire computational grid. If 
u s is the exact solution of the discretized equations, then 
ry = 0. In practice we obtain values for ry which are as 

small as roughly 10 -9 of the individual terms of f(u s )). 
This measure ensures convergence to the correct solution 
of the metric equations. 




60.0 



70.0 80.0 90.0 100.0 110.0 120.0 130.0 140.0 150.0 160.0 

t [ms] 

Fig. 2. The collapse phases used for specifying the met- 
ric resolution for regular collapse models (upper panel) 
and multiple bounce collapse models (lower panel). In a 
regular collapse, the metric is calculated on every nm^th 
time step in the pre-bounce phase (a). In the bounce phase 
(b), defined by p max > 0.1p nuc , the metric is calculated ev- 

(2) 

ery rim th time step. In the post-bounce phase (c), which 

(3) 

begins 5 ms later, the metric is calculated every rim th 
time step. In a multiple bounce model, the inter-bounce 
phases (a), during which the metric is calculated on ev- 
ery nm^th time step, are interspersed with bounce phases 

(2) 

(b), where the metric is calculated every n m th time step. 
Due to their lower average maximum densities, the density 
threshold is specified by /5 max > 0.02 / o nuc in these models. 



4.5. Metric extrapolation 

We have performed comparisons between the execution 
time for one metric solution (for different linear solver 
methods) and one hydrodynamic time step. The compar- 
ison shows that, particularly for grid sizes appropriate for 
our collapse simulations, the computation is completely 
dominated by the metric solver: One metric computation 
can be as time-consuming as about 100 hydrodynamic 
time steps! 

On the other hand, as long as the metric quantities 
do not change too rapidly with time, it is a fair approx- 
imation to solve for the metric not at every time step. 
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Table 3. Values for the metric refinement parameters / re f , 
pthrj Ai re f, and ro m for regular collapse models and for 
multiple bounce scenarios as used in our simulations (see 
text for more details). 



Collapse 


/ref 


Pthx 


At rcf 


n (1) 


n (2) 


(3) 






[10 14 g cm" 3 ] 


[ms] 








Regular 


0.1 


0.2 


5.0 


100 


10 


50 


Multiple 


0.02 


0.04 




100 


10 





We introduce the metric resolution parameter An m , such 
that mod (n, An m ) = 0. Then An m expresses the number 
of "hydrodynamic" time steps between two "metric" time 
slices where the metric is calculated. 

During a simulation, the metric resolution parame- 
ter should vary with time in a way which adapts to the 
specific phase of the collapse scenario. We split the core 
collapse evolution into three different phases: The pre- 
bounce phase lasts from the start of the evolution until 
the maximum density on the grid p max reaches a threshold 
value, which is a specified fraction of the nuclear density: 
/9 t h r = JrcfPnuc- During that phase the metric is calcu- 
lated every rim "hydrodynamic" time step. The subse- 
quent bounce phase lasts for a time At rc f, and the met- 

(2) 

ric is solved on every nfn th time slice. In the subsequent 
post-bounce phase, the metric resolution parameter is set 



to rim ■ 

In models showing multiple bounces, i.e. where the 
maximum density has multiple distinct peaks in time of- 
ten separated by several 10 ms, the evolution is split into 
a series of consecutive bounce phases interspersed with 
phases where the density is below the threshold marked 
by pthr- During the bounce phases the metric resolution is 

(2) 

given by ni, , and before the first bounce and in between 
bounces by r£ '. 

The collapse phases introduced above are sketched in 
Fig. |^, while the actual values for / re f and Ai re f used in our 
simulations are summarized for the two different collapse 
types in Table |[ 

In order to approximate the actual evolution of the 
metric on those "hydrodynamic" time slices where it is 
not calculated, the state of the metric at old "metric" time 
slices can be used to extrapolate the metric quantities for- 
ward in time. Fig. |§| shows the approximation of the actual 
metric calculated at every time slice by a metric which is 
calculated at every 10th time step, and is kept constant or 
is linearly or parabolically extrapolated in between. Tests 
with collapse models demonstrate that our choice of val- 
ues for the metric resolution parameters is appropriate to 
obtain the desired accuracy. Due to the superior accuracy 
of the 3rd order parabolic metric extrapolation scheme we 
use this method in all our simulations. 

5. Code tests 

We present now the tests we have used to calibrate our 
code. It is worth pointing out that the present code also 



1.03815 



1.03810 



1.03805 




30.19 30.20 30.21 30.22 30.23 30.24 

t [ms] 

Fig. 3. Approximation of the metric evolution by extrapo- 
lation. The evolution of the central value of the conformal 
factor <j) c during core bounce calculated at every time step 
(solid line) is very well approximated by a parabolic ex- 
trapolation of the metric which is calculated at every 10th 
time step (dashed line, +). If the extrapolation scheme is 
linear (dotted line, x) or constant (dashed-dotted lines, 
*), the approximation is less accurate. The symbols mark 
the instants when the metric is actually calculated. The 
model used for this plot is A1B3G5 in low grid resolution. 
For the usual resolution of 200 x 30 grid points the devia- 
tions between the different metric extrapolation schemes 
are hardly visible. 



includes the possibility of performing simulations of rota- 
tional core collapse employing Newtonian gravity. Many 
routines and numerical algorithms are common to both 
the relativistic and the Newtonian computer code, and by 
comparing with previous Newtonian simulations we have 
also checked the correct implementation of those routines. 
For this purpose we have run some of the models of the 
mprehensive sample of ^werger fc Miiller (1997 ) find- 



co 

ing excellent agreement for both the dynamics of the col- 
lapse and the gravitational waveforms. Therefore, when- 
ever we present Newtonian results below, it is assumed 
that those have been computed with our current code. 
We focus hence on checking our code in the relativistic 
regime. 

5.1. Relativistic shock tube tests 

Shock tubes can efficiently test the ability of a code to 
handle shock fronts. In a shock tube test, initial data are 
given for a fluid system which has two states of constant 
pressure, density and zero velocity, initially separated by a 
boundary. The time-evolution of these initial data yields a 
combination of constant states separated by shocks, con- 
tact discontinuities, or rarefaction waves. These tests al- 
low to compare the numerical results against the analytic 
solution flMartf fc Miiller 1994| ) in a straightforward way. 
We have simulated the two relativistic shock tube 



problems analyzed by Marti & Miiller (1996). The ini- 
tially constant states are given in Table |4j. Both cases 
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Table 5. Parameters of the rotating neutron star model 
used for the stability test. p c is the central density, 7 and 
K specify the EoS, r c /r p is the ratio of equatorial to po- 
lar radius of the star, Af rest / grav are the rest/gravitational 
masses, /? ro t is the rotation rate, fl/f^K is the ratio of the 
star's angular velocity to the Keplerian angular velocity 
and T is its rotation period. 



Pc 


7 


K 


r e /r p 


M rest 




p Iot si/Ok 


T 


[Pnuc] 




(cgs) 




[Af ] 


[M @ ] 


[%] [%] 


M 


3.952 


2.0 


1.456 x 10" 


0.70 


1.756 


1.627 


7.419 76.0 


1.0 



Fig. 4. Numerical (symbols) and analytic (solid line) pro- 
files of the (rescaled) density p, pressure P, and velocity 
v x for the shock tube problem 2. 

Table 4. Initial values of density p and pressure P for the 
left and right states of the relativistic shock tube prob- 
lems 1 and 2 of pVlarti fc Miiller (1996] ). 





Problem 1 
left state right state 


Problem 2 
left state right state 


p 
p 


10.0 1.0 
13.3 0.66 x 10" 6 


1.0 1.0 

10 3 10~ 2 



are computed in a flat one-dimensional Minkowski space- 
time using Cartesian coordinates in the interval x € [0, 1] 
spanned by 500 equidistant zones. The EoS is assumed to 
be that of an ideal fluid with adiabatic index 7 = 5/3, 
and th e initial discontinuity is placed at x — 0.5. The re- 
sults fc |r the most stringent problem 2 (a relativistic blast 
wave) are displayed in Fig. ||, after an evolution time of 
ifinai = 0.25. The profiles of the density p, the pressure 
P, and the velocity v x accurately reproduce the analytic 
solution. 

5.2. Rotating neutron stars 

An important test that any axisymmetric hydrodynamics 
code s hould pass is the computation of the evolution of 



fixed (Cowling approximation), or evolve both the hydro- 
dynamics and the spacetime in a fully coupled evolution. 
In both cases we observe that the neutron star remains in 
equilibrium to high accuracy for several 10 milliseconds, 
which corresponds to 10 rotation periods. In this test run, 
the metric equations are solved every 200th hydrodynamic 
time step, which corresponds to a time interval between 
two metric time slices of 0.06 ms. This interval is small 
compared to both the neutron star's rotation period and 
the frequencies of its strongest radial oscillations. 

In Fig. U we show the time evolution of the density 
p, the lapse a and the radial velocity v r — \J ^j xl v\ at 
a radius of 8.77 km in the equatorial plane for a fully 
coupled evolution. The small amplitude oscillations are 
triggered by numerical discretization errors. The fact that 
these oscillations are hardly damped reflects the low nu- 
merical viscosity of the HRSC scheme. By performing a 
Fourier transform of these evolved data one can obtain the 
frequen cies of the different mod es of pulsation to high ac- 
curacy flFant et al. 2000 , 2001a ). Comparisons perfo rmed 
with the resu lts in the Cowling approximation of Font 



et al. (2001a ) yield a very good agreement for the mode 
frequencies. 

As demonstrated in Fig. ||, the code can keep the equi- 
librium initial data to a high degree. The small secular 
drift observed in p and a is an artifact of the numerical 
scheme and has also been observed by Font et al. (2000). 
This drift can largely be minimized using a polytropic 
zero-temperature EoS, which is a fair assumption for 
studying small amplitude pulsations of polytropic stars. 



The effect of different EoS on the drift is discussed in Font 



an equ ilibrium model of a rotating stellar configuration. 
If the initial model is computed accurately, and if the evo- 
lution scheme is implemented properly, the matter and 
metric profiles of the evolved model should at most os- 
cillate slightly around their initial equilibrium states, i.e. 
the model should remain in equilibrium. In particular, for 
a rotating model, the rotation profile has to be preserved 
for as many rotation periods as possible. 

In order to enhance the relativistic effects we have used 
a rapidly and uniformly rotating neutron star model with 
a polytropic matter distribution (see Table ||) . This test 
allows for an independent check of the metric and the hy- 
drodynamics parts of the code. If desired, the code can 
solve only the fluid equations while keeping the metric 



et al. (2001b) 



In Fig. |6 we plot the radial profiles of the density p 
and the rotation velocity v v — \J 7 s3 V3 along the equator 
for the same neutron star model. Even after 5 rotation 
periods, the rotation profile is very close to its initial dis- 
tribution. Only at the stellar surface, the angular velocity 
slightly deviates from its initial shape due to the interac- 
tion with the atmosphere. As shown in Font et al. (2000| ) 
the use of high-order schemes such as PPM is essential in 
preserving the rotational profile. 

Since the neutron star model specified in Table |^ is 
rapidly rotating and thus deviates strongly from spherical 
symmetry, its evolution provides a convincing way to test 
the ability of the CFC to yield a good approximation of the 
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2.345 



2.340 




Fig. 5. Time evolution of hydrodynamic and metric quan- 
tities of the rotating neutron star model of Table || The 
density p (upper panel), the lapse function a (middle 
panel), and the radial velocity v r (lower panel) are dis- 
played at a radius r = 8.77 km in the equatorial plane 
(6 = tt/2). 



exact spacetime (see also Section 5.5). The initial model 
has been calculated using the exact metric, Eq. (^6|) , with- 
out assuming conformal flatness. The initial data arc then 
matched to the CFC spacetime. The fact that the evo- 
lution code maintains stability for many rotation periods 
already proves the validity of the CFC for rotating neutron 
star spacetimes. This observation is supported by Fig. [?] 
which shows how well the conformal factor <fp ¥c of the CF 
metric approximates the corresponding metric component 
4> ex of the exact metric. The excellent approximation of the 
azimuthal component of the exact shift vector f3 3 ex by the 
CFC shift vector component f3 3 is shown in the lower 
panels. Our result s are in good agreement with those of 



Cook ct al. (1996| ) 



5.3. Spherical core collapse 

The problem of relativistic core collapse for an ideal fluid 



so 




in spherical symmetry was first considered by May & 



White (1966). Their code used a Lagrangian formulation, 
where the coordinates label mass shells rather than spa- 



r [km] 

Fig. 6. Density profile (upper panel) and angular velocity 
profile (lower panel) of the rotating neutron star model of 
Table |5j. The equatorial density profile p c at t = 5.0 ms 
(dashed line), corresponding to 5 rotation periods is close 
to the initial profile at t — 0.0 ms (solid line). The same 
holds for the angular velocity profile v vc . The spike in 
the t = 5.0 ms profile of v^ c close to the neutron star's 
surface is a numerical artifact due to the artificial atmo- 
sphere. The height of the spike oscillates in time with a 
maximum amplitude as plotted here. Note that the neu- 
tron star model rotates at more than 30% of the speed of 
light at the equator. The dotted profile in the upper panel 
gives the initial density along the polar axis p p . Due to 
the rapid rotation, the polar radius of the neutron star is 
only 70% of its equatorial radius. The horizontal dotted 
line in the density plot marks nuclear matter density. 



tial positions and thus comove with the fluid. Artificial 
viscosity terms were included to damp out spurious nu- 
merical oscillations behind shock fronts. In this Section 
we compare our Eu lerian code against a Lagrangian fi- 
nite difference code ( Dimmelmeier 1998 ) for a spherically 
symmetric core collapse. 

For this test we have set up a nonrotating equilibrium 
star as initial data on a two-dimensional (r, #)-grid, with 
the same EoS and central density as in the (rotating) ini- 
tial models listed in Tab le fl} During the evolution, we do 
not use the hybrid EoS (|25|), but a simple ideal fluid EoS 
given by Eq. (p3|). 

The effective stiffening of the EoS at supranuclear mat- 
ter densities is modeled by a variable adiabatic index given 
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12.0 



20.0 



20.0 



Fig. 7. Quality of the CFC approximation for rapidly ro- 
tating neutron star initial data. The values for the exact 
and approximate conformal factor 0° x and (j)^ FC (upper 
top panel), and for the exact and approximate shift vec- 
tor component /3^ cx and (3^ CFC (lower top panel), both 
evaluated along the equatorial plane, agree very well. In 
the lower panels the corresponding relative deviations are 
plotted. For the conformal factor the relative deviations 
are less than 0.2%, while they stay below 4% for the shift 
vector. 



by the following relation ( van Riper 1979| ; Romero et al 
19960 



7 



7i 
7i 



-log — 

^ Pnm 



for p < p„ 
for p > p„ 



(56) 



with 71 = 1.320. Besides comparing our results with those 
of an independent spherically symmetric code, this test 
also allows us to assess the capability of the code to main- 
tain spherical symmetry throughout the collapse. 



6XJ 




69.0 



1000.0 



Fig. 8. Spherically symmetric core collapse: Comparison 
of results obtained with our Eulerian HRSC code (solid 
lines) against results from a Lagrangian finite difference 
code (dashed lines). The time evolution of the central den- 
sity p c (upper panel), here plotted at the time of bounce, 
agrees very well during all stages of the collapse. The dot- 
ted line in the density plot marks nuclear matter density. 
In the lower panel, the radial velocity profiles are plotted 
at t = 60 ms (a) and at t = 75 ms (b), respectively. During 
the infall phase before core bounce (a), the velocity pro- 
files match very well. When the shock front propagates 
outward (b), the agreement is less good due to the more 
dissipative character of the Lagrangian code. The inset 
shows the growth of unphysical oscillations in the veloc- 
ity profile obtained with the Lagrangian code for very low 
artificial viscosity at t — 75 ms. 



The differences in the ability of both codes to han- 
dle shocks are demonstrated in Fig. ^[ In the upper panel 
we plot the evolution of the central density around the 
time of bounce. We note that due to the relatively soft 
supranuclear EoS used for this run the core dives deeply 
into the potential well reaching densities as high as 5p nuc . 
Even during this epoch, where the highest densities are 
reached, the evolution computed with our code follows 
very accurately the corresponding one obtained with the 
May-White code. However, due to its superior effective 
radial resolution, the May- White code resolves the drop 
in the central density after the maximum peak better, and 
also produces a smoother density evolution. Up to the for- 
mation of the shock front, both codes yield equally good 
results for the radial infall velocity profiles. Such profiles 
are plotted in the lower panel of Fig. || at t = 60 ms (be- 
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fore bounce) and t = 75 ms (after bounce), respectively. 
With the HRSC scheme of our code the shock front is 
steeper and confined to fewer grid zones than with the 
artificial viscosity scheme implemented in the May- White 
code. Additionally, the velocity behind the shock is higher 
and the shock propagates faster, which is a result of much 
less numerical dissipation. We note that if the artificial 
viscosity is reduced in the May- White code to cure these 
negative effects, unphysical spurious oscillations grow be- 
hind the shock front (see inset in the lower panel of Fig. ^J). 

5.4. Integral quantities - Conservation of rest-mass 
and angular momentum 

In the 3 + 1 formalism of general relativity, the definitions 
for the (baryon) rest mass M rest , the proper mass M propor , 
the gravitational mass M grav , the angular momentum J rot 
and the rotational mass (energy) M rot are ( Friedman et al 
1986| ): 



M iest = - ^ 


1 pu^dV, 


(57) 


-^proper ~~ 


f p(l + e)u^dV, 


(58) 


-^grav 


f (2T flu -g^T)t"n v dV, 


(59) 


Jrot 


f T^n v dV, 


(60) 


M rot = - 


f rp n H n v S ^ U Ay 


(61) 



where V* and s M are a time-like and a space-like Killing 
vector, respectively. The rotation axis is perpendicular to 
s' 1 . For the CF metric in axisymmetry, where = e 3 
(i.e. the unit vector in the (^-direction), the above integral 
quantities have the following form: 



M rest = -2tt / r 2 sin e<j? pWdrdO, 



M erav = — 27T 



2tt J r 2 sm8<f) 6 p(l + e)Wdrd6, 

2 sin 6<j) 6 (a{2phW 2 + 2P - ph) 
2phW\/3 i ) drd6», 



J rot " " 



Mr, 



-2tt 
-2tt 



' sin 6<t> 6 phW 2 v 3 drdd, 



(62) 
(63) 

(64) 
(65) 



r 2 sin 0<f> 6 phW 2 Va ^ ^ drdO. (66) 



Note that for c = 1, the energies have the same units as 
the corresponding masses. The potential (mass) energy is 
given by M = M grav - M proper - M Iot . 

This formulation of the above integral quantities is an 
extension of the masses and the angular momentum de- 



fined in Komatsu et al. (1989a) for nonzero radial and an- 
gular velocity and shift vector components, applied to a 
conformally flat metric. It is obvious that the integrand of 
the rest mass and the angular momentum can be written 



as and y^S^, respectively. These are the first and 

fourth component of Eq. (|lj) and correspond to the con- 
tinuity equation and the angular momentum conservation 
equation, respectively. We have used these equations to 
check the conservation properties of our numerical code. 

First, we have performed tests where we have evolved 
the hydrodynamic state vector in a fixed spacetime met- 
ric. Additionally, the fluxes at the outer boundary are ex- 
plicitely set to zero and no artifical atmosphere is used. 
In this situation, we obtain conservation of the rest mass 
and the angular momentum of 



dM„ 



dJ T , 



M, 



Jrc 



rest 



l 



10 



-15 



10 



-15 



(67) 



(68) 



i.e. up to machine precision. Here M rest o and J ro to are 
the initial total rest mass and total angular momentum, 
respectively. In Newtonian gravity, the same level of con- 
servation is reached. 

In a dynamic core collapse simulation the above re- 
strictions have to be abandoned. The nonstationarity of 
the metric and (in the case of rotation) the artificial at- 
mosphere prevent a highly accurate conservation of the 
integral quantities. In Fig. || we plot the behavior of the 
mass and the angular momentum in the rotational core 
collapse model A3B2G4 (see Table 0). The dashed lines 
correspond to a Newtonian simulation and the solid ones 
are the relativistic counterparts. In Newtonian gravity the 
only source for mass and angular momentum "accretion" 
is the artificial atmosphere. This becomes visible as a 
small monotonous increase in M rcst (middle panel) and 
Jmt (bottom panel). Until bounce, indicated by the ver- 
tical dotted line, the mass increase in relativistic gravity 
exhibits the same behavior. However, due to the nonzero 
source term in the relativistic angular momentum equa- 
tion, dJ TO t is larger than its Newtonian counterpart al- 
ready during the infall phase. 

Around the time of bounce another effect which pro- 
hibits exact numerical conservation of total rest mass and 
angular momentum becomes important. The conservation 
equations ( [l8"|) are formulated for y/jU. But in our numer- 
ical scheme we evolve only the state vector U and, for nu- 
merical reasons, we treat the geometric volume term ^7 
as a source term. Therefore, for any strong temporal vari- 
ation of the metric, e.g. during bounce, dM les t and dJ IQ t 
vary considerably, too. This is clearly visible in Fig. ^[ 
When the central density has settled down to an almost 
constant value, and the core has reached a new equilib- 
rium configuration, the rates of increase of M res t and J rot 
are again similar in relativistic and in Newtonian gravity. 

The effects of the artificial atmosphere on the conserva- 
tion of total mass and angular momentum can be infered 
from Fig. |l0|. In the highly nonspherical rotational core 
collapse model A4B2G2, the shock propagates outward at 
different speeds in the polar and the equatorial direction. 
When the shock reaches the surface of the star at the 
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60.0 



Fig. 9. Conservation of total rest mass and angular mo- 
mentum in the relativistic (solid lines) and Newtonian 
(dashed lines) simulation of model A3B2G4. Until the 
time of bounce at % w 39 ms (indicated by the verti- 
cal dotted line) the relative changes in the total rest mass 
M^cst (middle panel) and angular momentum J rot (lower 
panel) due to interaction of the stellar interior with the 
artifical atmosphere are small compared to the respective 
initial values M res to and J ro to- Nonzero source terms in 
the evolution equation for S3 drive additional nonconser- 
vation of the angular momentum in relativistic gravity 
(indicated by the arrow). As the central density p c (upper 
panel) and the metric volume factor ^7 vary most during 
bounce, a strong oscillatory variation of total mass and 
angular momentum can be observed in this phase (asso- 
ciated with the pulsations observed in the proto-neutron 
star) . The horizontal dotted line in the density plot marks 
nuclear matter density. 




100.0 



Fig. 10. Influence of the artificial atmosphere on the con- 
servation of total rest mass and angular momentum for 
model A4B2G2. Before bounce the change of total mass 
and angular momentum due to interaction of the stellar 
interior with the atmosphere is small. Around the time 
of bounce at « 68.5 ms (a), the maximum density 
Pmax (upper panel), which is off-center in this toroidal 
model, and the volume element ^/y oscillate strongly. At 
t k 85 ms (b), the shock reaches the stellar surface, first 
at the pole at i?b P ~ 1 200 km. This leads to mass loss set- 
ting in around that time. At t « 92.5 ms (c), the slower 
shock along the equatorial plane breaks through the stel- 
lar surface, which is located at Rh c ~ 1320 km > i?b p . 
This results in a loss of total angular momentum, and also 
increases the mass loss. 



the variation of dM res % and dJ rot due to rapid changes of 
the volume element x/7 around peak maximum density. 



pole (t ss 85 ms), a considerable amount of mass, but not 
much angular momentum is lost to the atmosphere. The 
equator is reached by the shock at a later time, since the 
shock propagates slower in this direction and the equato- 
rial radius is larger. Around that time at t w 92.5 ms the 
mass loss curve steepens and significantly more angular 
momentum is lost, as the relative contributions of fluid 
elements to the total angular momentum is higher near 
the equatorial plane than at the pole. Fig. also shows 



5.5. Quality of the conformally flat condition during 
rotational core collapse 



As mentioned in Section 5.2, the accuracy of the CFC ap- 



proximation has been considered by Cook et al. (1996), 



who find remarkably good results for rapidly rotating rel- 
ativistic stars in equilibrium, typical errors for different 
variables being smaller than 5%. However, the quality of 
the CFC approximation degrades in the case of extremely 
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relativistic nonspherical configurations, e.g. for rigidly ro- 
tating infinitcsirnally thin disks of dust, which have been 
investigated by Kley & Schafer (1999). These authors also 
point out that general relativity and the truncated CFC 
approach are identical up to the first post-Newtonian or- 
der. Although the above comparisons have been made for 
stationary configurations, we want to emphasize that the 
CFC approximation does not make any explicit assump- 
tions about the spacetime being static or dynamic. Due 
to the elliptic character of the metric equations ([?] ||), 
the spacetime metric is solely determined by the instanta- 
neous hydrodynamic state. Therefore the deviation of the 
CF metric from the exact spacetime metric at a particu- 
lar time will only be due to the current deviation of the 
matter distribution from spherical symmetry. 

In this Section we analyze the validity of the CFC in 
simulations of rotational core collapse. There are several 
ways to quantify the violations of the exact metric equa- 
tions. Firstly, one could use the CF metric obtained by 
solving Eqs. to compute the Einstein tensor G^ u and 
to check the Einstein equations = SitT^. However, 
this is a difficult approach due to the huge number of terms 
involved. 

On the other hand, the ADM metric equations 
are equivalent to the Einstein equations. As we only use 
a subset of the ADM equations to calculate the CF met- 
ric, we can use the remaining ones to test the quality of 
the CFC approximation. Therefore, we re-insert the CF 
metric components into the individual ADM metric equa- 
tions and compare the left and right hand sides of the 
respective equations. We note that such a procedure has 



also been used by Gourgoulhon et al. (2002) in validat- 



ing initial data for binary black holes. Comparisons with 
Damour's analytic post-Newtonian expansion show satis- 
factory agj^e^mcn^hithis case up to third post-Newtonian 



order ( Gourgoulhon 2001 

If the inserted metric is an exact solution of the 
Einstein equations, the numerical inequality between the 
left and right hand sides of the ADM equations has to 
converge to zero with increasing resolution. For an ap- 
proximate metric like the CF metric there will remain a 
nonzero residual 



lhs q 



rhs q 



1 



> 



(69) 



in the ADM equations for at least some metric quantities 
q irrespective of numerical resolution. Here lhs q and rhs q 
stand for the left and right hand side of the corresponding 
ADM metric equation, respectively. 

We want to point out that the metric equations (0-[|) 
are analytically equivalent to the ADM constraint equa- 
tions. Any solution of the CF metric equations will also 
be a solution of the ADM constraint equations, at least 
up to the precision of the numerical scheme. 

The evolution equations for the three-metric jij, 
Eq. (Q), are also identically fulfilled for the off-diagonal 
elements of 7^ , because for a CF metric they transform 
into the definition of the extrinsic curvature, Eq. (113). As 
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Fig. 11. Numerical equivalence of the left hand side 
(lhs 4> — dt4>) and right hand side (rhs <j> = |0V/c/3 fe ) of the 
evolution equation for the conformal factor (f>. The resid- 
ual r§ of the evolution equation for 0, (plotted along the 
equatorial radius at t « 41 ms) is small everywhere. Note 
that the pole at large r is due to a vanishing denominator 
in Eq. ©. 



the three-metric 7y is conformally flat, of the six evolution 
equations (|^) there remains only the evolution equation for 
the conformal factor 0, Eq. ( pd| ) which must be monitored. 
We have found that the residual of this equation, 



lhs 4> 




rhs cf> 





Or 



1 



(70) 



is very close to zero in various test situations including ro- 
tational core collapse. For model A3B4G3 the radial pro- 
file of is plotted in Fig. ^ shortly before maximum 
density is reached. Except where a pole develops due to a 
zero denominator in Eq. ( |70|) , is everywhere less than 
0.2%, particularly in the dense interior. Hence, we infer 
that the ADM evolution equations for the diagonal ele- 
ments of 7ij are solved to high degree of accuracy by the 
CF metric as well. 

The ADM evolution equations for the extrinsic cur- 
vature provide another way of checking the accuracy of 
the CFC approximation. As a consequence of the max- 
imal slicing condition, K — 0, the sum of the left and 
right hand sides of the equation for the diagonal elements 
of the extrinsic curvature are identical up to numerical 
discretization errors at all grid points in any evolution. 
Furthermore, in spherical symmetry, the left and right 
hand sides of all six evolution equations for Kij agree, be- 
cause the CF metric is an exact solution of the spacetime 
in this case. 

A violation of the evolution equations for the diagonal 
elements of Kij, Eq. (^|), occurs in all highly nonspherical 
and relativistic spacetimes, e.g. in the core collapse model 
A3B4G3 (Fig. |l2] ). The agreement between both sides of 
the equations is good at low densities at the beginning 
of the infall phase. However, when high densities are en- 
countered in the center during and after bounce the mis- 
alignment becomes as prominent as for the rotating equi- 



librium neutron star model of Section 5.2 (bottom panel 
of Fig. 12). On the other hand, the evolution equations 
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Fig. 12. Effect of the CFC approximation on the ADM 
evolution equation for if 13 (upper panel) and ifn (lower 
panel) in the rotational core collapse model A3B4G3. In 
the case of if 13, the equatorial radial profiles of the left 
hand side (solid lines) and of the right hand side (dashed 
lines) of the ADM evolution equation agree very well down 
to small radii where relativistic effects are most important. 
For if 11 the terms — ViVj-a and aRij are nonzero. This 
causes the non-negligible misalignment between both sides 
of that equation (bottom panel). The comparison is done 
at t = 41 ms, which is about 1 ms before bounce. 

for off-diagonal components are solved to high accuracy. 
This is demonstrated by he very good agreement of the 
left and right hand sides of the equation for if 13 even in 
the center of the core (upper panel of Fig. |l^) . 

These observations raise the question why the left and 
right hand sides of some of the evolution equations for the 
extrinsic curvature disagree considerably for nonspherical 
matter distributions. The reason becomes obvious when 
analyzing the structure of those equations. In comprehen- 
sive tests computing many dynamic collapse models we 
have found that the terms in Eq. (^| containing extrinsic 
containing extrinsic curvature components arc negligibly 
small compared to the other three terms — ViV^a, aRij 
and — 8na(Sij — ^lij{S — ps)) everywhere on the grid and 
during all epochs. The relevant terms depend only on the 
metric quantities a (which is fixed by a gauge condition) 
and cj), but not on the shift vector components f} 1 . This 



points towards the inability of the CF three-metric 7^ , 
which has only one nonzero component <fi, to exactly rep- 
resent the physical spacetime. On the other hand, we note 
that the quality of the shift vector components, which are 
calculated directly according to the ADM momentum con- 
straint equations (|5|) (which are equivalent to the metric 
equation (||)), is not measured by Eq. (||). 

Therefore, the misalignment of the evolution equation 
for the extrinsic curvature is not a result of strong devi- 
ations of the shift vector components, but is instead due 
to the restriction imposed onto the three-metric 7^ by 
the CFC approximation. This conjecture is additionally 
supported by the fact that, as demonstrated in Fig. [?] in 
Section 5/2, even for high density, rapidly rotating neutron 
stars in equilibrium the value of the shift vector compo- 
nent (3 3 calculated by the CFC equation @ remains close 
to the exact value. 

This also explains the perfect matching of the left and 
right hand sides of Eq. (J3J) for if 13 even for rapidly rotat- 
ing and collapsing models. For this component the prob- 
lematic terms — V^V^a and aRij vanish. However, for if 12 
and if 23 the deviations are more pronounced, since the po- 
lar velocity vi is very close to zero in our collapse models. 
As a result the left and right hand sides of these evolution 
equations are dominated by truncation errors. 



6. Conclusions 

We have presented a detailed description of a new axisym- 
metric general relativistic code for rotational core collapse 
which solves the coupled system of Einstein equations and 
fluid equations. Those equations are integrated adopting 
the ADM 3+1 formalism to foliate the spacetime into a 
set of non- intersecting spacelike hypersurfaces. In order 
to simplify the metric equations we approximate the gen- 
eral metric usi ng the conformal fla tness condition (CFC) 
introduced by Wilson ct al. (199(: ), whereby the spatial 
components of the metric are set equal to the flat three- 
metric times a conformal factor which depends on the co- 
ordinates. The CFC implies a considerable simplification 
of the ADM 3+1 metric equations which reduce to a set 
of five coupled non-linear elliptic equations for the metric 
components. We have carefully analyzed the applicability 
and quality of the CFC, showing that this approximation 
is appropriate for simulations of rotational core collapse. 

We have also demonstrated that our general relativistic 
hydrodynamics solver accurately passes a set of difficult 
test problems including special relativistic shock tubes, 
the preservation of the rotation profile and of the equi- 
librium of rapidly and differentially rotating relativistic 
polytropes, and the spherical collapse of a relativistic poly- 
trope. The solver, which is based on a high-resolution 
shock capturing scheme with an approximate Riemann 
solver, also conserves accurately rest-mass and angular 
momentum in dynamic spacetimes. 

The code described in this work was developed with 
the aim to simulate the general relativistic axisymmetric 
collapse of rotating stellar cores. To this end we have com- 
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puted and presented here a large set of relativistic rotating 
precollapse equilibrium models with different amounts and 
distributions of angular momentum. Using these initial 
models we have performed a series of collapse simulations 
and calculated the resulting gravitational wave signal. The 
results of this first relativistic study of rotational core col- 
lapse and of its gravitational wave signature are discussed 
in an accompanying paper (Dimmelmeier et al. 2002). 



Acknowledgements. We are grateful to Y. Eriguchi for pro- 
viding us with a code for constructing the intial models. It 
is a pleasure to acknowledge J. M- Ibanez and G. Schafer 
for very useful comments and suggestions. We also want to 
thank M. Rampp for insightful discussions about the sweep- 
ing method for sparse linear problems. J.A.F. acknowledges 
support from a Marie Curie fellowship by the European 
Union (MCFI-2001-00032) and from the Spanish Ministerio 
de Ciencia y Tecnologfa (AYA2001-3490-C02-01). All compu- 
tations were performed on the NEC SX-5/3C supercomputer 
at the Rechenzentrum Garching. 

References 

Arnowitt, R., Dcser, S., & Misner, C. W. 1962, in 

Gravitation: An introduction to current research, ed. 

L. Witten (New York, U. S. A.: Wiley), 227-265 
Banyuls, F., Font, J. A., Ibanez, J. M., Marti, J. M., & 

Miralles, J. A. 1997, Astrophys. J., 476, 221 
Colella, P. & Woodward, P. R. 1984, J. Comput. Phys., 

54, 174 

Cook, G. B., Shapiro, S. L., & Teukolsky, S. A. 1996, Phys. 

Rev. D, 53, 5533 
Dimmelmeier, H. 1998, Diploma thesis, Universitat 

Regensburg, Regensburg, Germany 
Dimmelmeier, H., Font, J. A., & Miiller, E. 2001, 

Astrophys. J. Lett., 560, L163 
— . 2002, Astron. Astrophys., submitted 
Donat, R. & Marquina, A. 1996, J. Comput. Phys., 125, 

42 

Flanagan, E. E. 1999, Phys. Rev. Lett., 82, 1354 

Font, J. A., Dimmelmeier, H., Gupta, A., & Stergioulas, 

N. 2001a, Mon. Not. R. Astron. Soc, 325, 1463 
Font, J. A., Goodale, T., Iyer, S., et al. 2001b, Phys. 

Rev. D, in press 
Font, J. A., Stergioulas, N., & Kokkotas, K. D. 2000, Mon. 

Not. R. Astron. Soc, 313, 678 
Friedman, J. L., Ipser, J. R., & Parker, L. 1986, 

Astrophys. J., 304, 115 
Gourgoulhon, E. 2001, private communication 
Gourgoulhon, E., Grandclement, P., & Bonazzola, S. 2002, 

Phys. Rev. D, 65, 044020 
Janka, H.-T., Zwerger, T., & Monchmeyer, R. 1993, 

Astron. Astrophys., 268, 360 
Kley, W. & Schafer, G. 1999, Phys. Rev. D, 60, 027501 
Komatsu, FL, Eriguchi, Y., & Hachisu, I. 1989a, Mon. Not. 

R. Astron. Soc, 237, 355 
— . 1989b, Mon. Not. R. Astron. Soc, 239, 153 
Marti, J. M., Ibanez, J. M., & Miralles, J. A. 1991, Phys. 

Rev. D, 43, 3794 



Marti, J. M. & Miiller, E. 1994, J. Fluid Mech., 258, 317 
— . 1996, J. Comput. Phys., 123, 1 

Mathews, G. J. & Wilson, J. R. 2000, Phys. Rev. D, 61, 
127304 

May, M. M. & White, R. H. 1966, Phys. Rev., 141, 1232 
Monchmeyer, R., Schafer, G., Miiller, E., & Kates, R. E. 

1991, Astron. Astrophys., 246, 417 
Miiller, E. 1982, Astron. Astrophys., 114, 53 
Miiller, E. 1998, in Computational methods for astro- 
physical fluid flow. Saas-Fee Advanced Course 27, ed. 
O. Steiner & A. Gautschy (Berlin, Germany: Springer), 
343-494 

Potter, D. 1973, Computational physics (Chichester, 
U. K.: Wiley) 

Pradier, T., Arnaud, N., Bizouard, M.-A., et al. 2001, 

Phys. Rev. D, 63, 042002 
Romero, J. V., Ibanez, J. M., Marti, J. M., & Miralles, 

J. A. 1996, Astrophys. J., 462, 839 
Stergioulas, N. & Friedman, J. L. 1995, Astrophys. J., 444, 

306 

Toro, E. F. 1997, Riemann solvers and numerical methods 

for fluid dynamics - a practical introduction (Berlin, 

Germany: Springer) 
van Riper, K. A. 1979, Astrophys. J., 232, 558 
Wilson, J. R., Mathews, G. J., & Marronetti, P. 1996, 

Phys. Rev. D, 54, 1317 
Zwerger, T. 1995, PhD thesis, Technische Universitat 

Miinchen, Miinchen, Germany 
Zwerger, T. & Miiller, E. 1997, Astron. Astrophys., 320, 

209 



Appendix A: Recovery of the primitive variables 

For an ideal gas EoS a suitable method to recover the 
primitive quantities from the state vector was introduced 



by Marti et al. (1991). It has been described in detail in 
Marti fc Miiller (1996 ). We have extended this method to 
be adequate for our hybrid EoS (p5j). 

We calculate a function f(P*) = P(p*,e*) - P* , where 
p* and e* depend on the conserved quantities and P* only 
(in all relations involving the pressure we replace P by 
P*): 



f(P*) = K (Dy/xY ( 1 - 



+( 7th - 1)X 



7th " 1 
7-1 

t + D ( 1 - 



P* 1 

X 



(7th ~ 1)(7 - 7i ) 
(7i - !)(72 - 1) 



KpZ^dVx-p*, 



(A.l) 



where we use the auxiliary quantity X = 1 — S 2 /(t + 
D + P*) 2 , and S — SiS 1 . The new pressure is iteratively 
computed using the Newton-Raphson method: 



P* 



(A.2) 
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where the derivative of f(P*) is given by 

d/(P * } ir(Wx)V( 



H. Dimmelmeier et al.: Relativistic rotational core collapse. I 

4ii 3 iA 



1 _ 7th -1 
7-1 



dP* 



+ 



X(t + D + P*) 3 
2(7th-l)S 2 [ T + ^(l--U+P*(l-^) 



+(7th - l)X 



(t + D + P*) 
DS 2 



X 3 l 2 (r + D + P*) 3 
1 2P*S 2 



f [U ) = u H 1 r 



X X 2 (r + P + P*) 3 
(7th -1) (j-^Kpl^DS 



(7i-l)(72-l)V^(r + -D + P*) 3 



1. (A.3) 



During each iteration, the primitive quantities and the 
Lorentz factor, which are needed in the EoS, are updated 
using the following relations: 

* (A.4) 



W* 



T + D + P* 

1 

V 1 - v * v " ' 

D 

t + D(1- W*) + P*(l - W* 2 ) 
DW* ' 



(A.5) 

(A.6) 
(A.7) 



The iteration is continued until convergence is reached. 



Appendix B: Explicit form of the metric equations 



Using the notation of Eqs. ( p2|) and (|53|), the metric equa- 
tions read (with uj = dfU, uj g — dfd g u): 

tit fc A i , 2u i , u ]ee , cote l 
f (u )= u, rr + + -hr + —o-u d 



+ 2TT(u 1 ) 5 (phW 2 - P) 

ju 1 ) 7 
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+4«) 2 + -Sr'iu^y + 3 sin 2 9{u b fi f 

4:U 3 U 4 , 
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r 2 r r u 1 



u , 



1 

it 2 

\ {< 

cot 9u 4 
2u 3 



7u 2 v} s 
7uV r \ 



u 



u 



2ut - u , 



\r6 



2u 3 
r 



2u 3 



cot 9u 4 r 



4ut u 4 



0. 

cot 6 



4 (u fc ) = U 4 „ + — 2- + -^j- 

,4 1C^,2 



+ 



(l-C0t6»)M 4 leTTli 2 ^ 



r*2/i 1 1 



1 

« 2 

2 
S^ 2 



w 2 a - 



lu 2 u\ 
lu 2 u 4 B 



-u 3 r + 2u 4 g 



H cot 9u 

r 



1 



2u 3 



U 



sm 



4ul 



3r 
= 0, 

3cot0 



. I U ,r6 + W ,l 



r r 
16nu 2 S 3 



1 

r 2 sin 2 0m 1 u 2 
1 / „ 7uV, 



2 7u 2 u 4 r 



= 0. 



0. 



(B.l) 



?7/ 2 



If 



00 



fc\ 2 i *>T | ,oo 

/(■")= <rr + + -3- 



cot6> 



2;t(ii 1 ) 4 ii 2 (p/i(3V^ 2 - 2) + 5P) 



?( u ) 



3r 2 sin 2 6>(u 5 J 2 + 



3(i 



48(u 2 ) 2 

+4(u 3 r ) 2 + 3r 2 (7i 4 r ) 2 + 3 sin 2 9{u 5 e ) 2 



4(n 3 ) 2 
8it 3 it 3 



(B.2) 



(B.3) 



(B.4) 



(B.5) 



